
// separation of variables method
// With resampling....

template<class Resample>
class Orthantst_PF{
	public:
		Orthantst_PF(double *L, double *b, int p, int nu,Resample *R)
		{
			_L=L;
			_b=b;
			_p=p;
			_st=new boost::math::students_t_distribution<>(nu);
			_nu=nu;
			_R=R;
		}
		~Orthantst_PF()
		{
			delete _st;
		}
		void operator()(int M,double *res, double *eta)
		{
   			mat u(M,_p);
   			mat y(M,_p);
			double pa,pb,tpz,arg,tp1,ru;
			double sum2=0;
	   		double thres=*eta*M;
   			boost::random::uniform_01<> U;
   			RandomG::Random<boost::random::uniform_01<> > *_runif=new RandomG::Random<boost::random::uniform_01<> >(U);
   			double* prod=new double[M];
			for(int i=0;i<M;i++)
   			{
	  			prod[i]=1;
			}


			for(int j=0;j<_p;j++)
			{

				double ess=Ess(prod,&M);
				//    cout << "ess " << ess <<"\n";
				if(ess<thres){	
					 cout << "ess " << ess <<" ";
					*res+=mean2(prod,&M);
				//	cout << mean2(prod,n); 
					fResample(prod,&M,&y,&u);
				}
				for(int i=0;i<M;i++)
				{
					sum2=0;
					tp1=sqrt((as_scalar(sum(y(i,span(0,j))))+_nu)/(_nu+j));	
					for(int k=0;k<j;k++) sum2+=_L[k*(_p)+j]*u(i,k);
					tpz=(_b[j]-sum2)/(_L[j*(_p)+j]*tp1);
					pa=cdf(*_st,tpz);
					pb=1;
    					prod[i]*=pb-pa;
				//	cout << prod[i] << "\n";
					ru = (*_runif)();
		 			arg=ru*pb+(1.-ru)*pa;
					if(arg > .999999999) arg=.999999999;
					if(arg < .0000000001) arg=.0000000001;
					u(i,j) = quantile(*_st,arg);
					if(j+1<_p)
					{
						y(i,j+1)=u(i,j)*sqrt((as_scalar(sum(y(i,span(0,j))))+_nu)/(_nu+j));
					}
				}				
			}	
		  *res=prod[0];
		  for(int k=1;k<M;k++)
		  {
			  if(exp(prod[k])!=0)
			  {
				   *res=log_add(prod[k],*res);
			  }
		  }
		  *res=*res-log(M);
		   delete _runif;
		   delete[] prod;


		}

		inline void Normalize(double *w,int *n)
		{
			for(int i=0;i<*n;i++)
			{
				w[i]=1;
			}
		}
		inline double mean2(double *w, int *n){
			double sum=log(w[0]);
			for(int i=1;i<*n;i++){
				if(w[i]!=0)
				{
					sum=log_add(log(w[i]),sum);
				}
			}
			return sum-log(*n);
		}
		inline void fResample(double *w, int *n,mat *z, mat *U)
		{
			mat v=growingvect(*n);
			(*_R)(n,w,v);
			Arangemat(*z,v);
			Arangemat(*U,v);
			Normalize(w,n);
		}
	
		inline double Ess(double *w, int *n)
		{
			double sum=0;
			for(int i=0;i<*n;i++)
			{
				sum+=w[i];
			}
			double sum2=0;
			for(int i=0;i<*n;i++)
			{
				sum2+=pow(w[i]/sum,2);
			}
			return (double)1/sum2;

		}
		

	private:
		double *_L;
		double *_b;
		int _nu;
		Resample *_R;
		boost::math::students_t_distribution<> *_st;
		int _p;
};
